function interpolProf_f(np, x, y, no, xo) result(yo)

    use iasing_osse

    implicit none
    ! Arguments
    integer(JPI4), intent(in) :: np, no
    real(JPR4), dimension(np), intent(in) :: x, y
    real(JPR4), dimension(no), intent(in) :: xo
    real(JPR4), dimension(no) :: yo

    ! internal variables
    integer(JPI4)   :: i, jp, kmin
    logical         :: lreverse 
    real(JPR4), dimension(no) :: xaux
    real(JPR4) :: a, b
    real(JPR4) :: w1, w2

    include "dichotom_f.interface"
    !print*

    ! Check correct order in input x arrays
    if (xo(1)>xo(2)) then
        !print*, "Inverting profile"
        xaux = xo(no:1:-1)
        lreverse = .true.
    else
        lreverse = .false.
        xaux = xo
    endif
    
    ! Begin interpolation
    do i = 1, no
        jp = dichotom_f(np, x, xo(i))
        !if (i==no) then
        !    print*, "jp: ", jp
        !    print*, "xo(no)", xo(i)
        !endif
        if (jp.gt.0 .and. jp.lt.np) then
            ! Linear interpolation
            w1 = (x(jp+1)-xo(i)) / (x(jp+1)-x(jp))
            w2 = (xo(i)-x(jp))   / (x(jp+1)-x(jp))
            yo(i) = y(jp)*w1 + y(jp+1)*w2
        else if (jp==np) then
            yo(i) = y(np)
        else if (jp.eq.-1) then
            yo(i) = y(1)
        else if (jp.eq.-2) then
            ! Constant extrapolation below the lower level of the reference
            ! profile
            yo(i) = y(np)
        end if 
    enddo
    
end function
